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£h Abstract 

o; 

in the framework of mapped pseudospectral methods, we use a polynomial-type mapping function in order to describe accurately 
^2 the dynamics of systems developing small size structures. Using error criteria related to the spectral interpolation error, the 
C_) polynomial-type mapping is compared against previously proposed mappings for the study of collapse and shock wave phenomena. 

" ^ As a physical application, we study the dynamics of two coupled beams, described by coupled nonlinear Schrodinger equations 
^►^and modeling beam propagation in an atomic coherent media, whose spatial sizes differ up to several orders of magnitude. It is 

^ C~| demonstrated, also by numerical simulations, that the accuracy properties of the polynomial-type mapping outperform in orders 
f**> magnitude the ones of the other studied mapping functions. 

Key words: mapping function, Chebyshev approximation, pseudospectral methods, partial differential equations, nonlinear Schrodinger 
^ equation 
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1. Introduction 

The numerical simulation of physical systems which may 
develop multiple scale phenomena, like damage fracture, 
tumor growth, transport and flow in heterogeneous media, 
propagation of (non)linear waves, has to be handled with 
care in order to properly reproduce all their physical fea- 
tures. In such situations the size of the spatial grid and 
the time advancing step may become critical issues for cap- 
turing the dynamics of this type of systems. The numer- 
ical difficulties related to a naive increasing of the num- 
ber of discretization points can be overcome by using more 
sophisticated techniques, for instance, domain decomposi- 
tion [23], multi-scale finite element method (FEMs) [7] or 
transformations through changes of variables [5] . Domain 
decomposition split the original domain into smaller subdo- 
mains which are independently discretized but still linked 
together by their boundary conditions, which have to en- 
sure a sufficiently smooth solution across the non-matching 
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grids of the different subdomains. Multi-scale FEMs take 
advantage of the construction of a specific set of basis func- 
tions according to the spatial size of each element of the 
mesh. There is in fact a broad class of FEMs dedicated to 
the analysis of multiple scale phenomena, each method be- 
ing designed to address a specific issue, for example, one 
can capture the large scale behavior of the solution with- 
out resolving all the small scale features [16]. On the other 
hand, domain transformation methods (or mapping func- 
tions) make use of bijective applications to map the points 
of the physical domain into a computational domain where 
the function to be discretized is to show a much smoother 
behavior. 

The use of spectral methods has become popular in the 
last decades for the numerical solution of partial differen- 
tial equations (PDEs) with smooth behavior due to their 
increased accuracy when compared to finite-differences or 
finite-elements stencils with the same degree of freedoms. 
This happens because the rate of convergence of spectral 
approximations depends only on the smoothness of the so- 
lution, a property known in the literature as "spectral accu- 
racy" . On the contrary, the numerical convergence of finite- 
differences or FEMs is proportional to some fixed negative 
power of N, the number of grid points being used. 



For problems with a less smoother behavior, such as those 
exhibiting rapidly varying solutions, there is a great deal 
of computational evidence that appropriately chosen map- 
ping functions can significantly enhance the accuracy of 
pseudospectral applications in thse situations, thus avoid- 
ing the use of fine grids and their associated spurious conse- 
quences. Examples include mappings to enhance the accu- 
racy of approximations to shock like functions [1-3, 11, 18, 
25] , approximation of boundary layer flows in Navier-Stokes 
calculations [8], multidomain simulation of the Maxwell's 
equations [14], or cardiac tissue simulations [28]. There is 
also considerable computational evidence that the changes 
in the differential operator introduced by the mapping do 
not negatively affect the conditioning of the matrices ob- 
tained from the pseudospectral approximation [1-3, 10, 11]. 

In this work we use a two-parameter polynomial-type 
mapping function in order to simulate the propagation of 
two coupled electromagnetic beams of transverse widths as 
disparate as up to three orders of magnitude. The parame- 
ters of the mapping function are adjusted in order to min- 
imize functionals related to the spectral interpolation er- 
ror. The polynomial mapping is compared against two pre- 
viously proposed mappings for shock-like fronts and wave 
collapse phenomena [4, 27] . 

The paper is organized as follows. In Section[2]we give 
a brief description of the underlying physical system. In 
Section [3] the polynomial mapping together with the other 
mappings are compared using error criteria, and the differ- 
ences between them are pointed out. In Section [4] the nu- 
merical scheme is presented and simulations of the phys- 
ical system arc performed using each mapping. Finally, 
Section [5] briefly summarizes our main conclusions. 

2. Physical system 

Atomic coherent media were brought into the focus of 
the scientific community with the theoretical proposal 
and experimental demonstration of electromagnetic in- 
duced transparency (EIT) [13]. EIT phenomena consists 
in rendering transparent a rather opaque media by means 
of an external electromagnetic field, and it is the result 
of destructive interference between two transition paths 
having the same final state [13]. The atomic coherent me- 
dia exhibits far more physical phenomena [24], like lasing 
without inversion, huge enhancement of refractive index, 
or negative refractive index [17]. 

The atomic coherent media of our interest is modeled by 
a noninteracting atomic gas possesing the four-level energy 
diagram shown in Fig. [T^,. The atom-fields interaction in- 
cludes the following parameters: relaxation rates 713, 723, 
724, the decoherence rate 712 between levels |1) and |2), 
the amplitudes of electromagnetic fields ^13, SI23, and 
the detunings A 13, A23, A24 of the field frequency with 
respect to the energy levels of the atomic media. A more 
detailed presentation of our four-level system can be found 
in Ref. [15] and the references therein. 
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Fig. 1. Diagram of the energy levels of the atomic system (a), and the 
dependence of nonlinear susceptibilities experienced by the probe, 
Xp, (b) and coupling, Xo (c) fields. 



Assuming an instantaneous response of the atomic me- 
dia to the electromagnetic fields, the beams propagation is 
modeled by a system of two coupled, two-dimensional non- 
linear Schrodinger (NLS) equations 

i^ = -Afi p - Xp (|tt p | 2 ,|tt c | 2 )r! p (la) 
i^ = -AO c - Xc (|ft p | 2 ,|tt c | 2 R, (lb) 

where fl p and f2 c are respectively known as the probe and 
coupling (control) fields, and \v an d Xc are the nonlinear 
susceptibilities of the atomic media experienced by these 
probe and coupling fields, respectively. In general, these 
susceptibilities exhibit both real and imaginary parts. For 
simplicity, in the present work we neglect the imaginary 
parts, which are actually associated with the fields absorp- 
tion. The susceptibilities can then be written in analytical 
form as the quotient of two bilinear forms of arguments 
I Op 1 2 and |f2 c | 2 , and are similar in structure to those de- 
rived in Ref. [26]: 

Xp ' c EijhiWW* n T p -B.U c ' u 

where n£ c = [1 |f! p , c | 2 \^c\ 4 l^ P ,c| 6 ■ ■ ■ |^ c | 2m - c ], with 
m p = 6 and m c — 5, are vectors, and A( P:C ) = {af- c ^}, 
B = {bij} are (m p + l) x (m c + l) matrices. The coefficients 
of these matrices are sensitive to the values of the fields de- 
tunings A 13, A23 and A24. For our particular configuration 
of fields detunings (712 = 1CT 8 7, 713 = 723 = O.67, 724 = 
1.257 and A13 = A23 = A24 = 57, where 7 = 30MHz is a 
normalization constant) matrices A p , A c , and B are given 
below. This configuration of detunings was motivated by 
the cubic-quintic-like model of the NLS equation, which can 
display liquid light behavior [19,20]. In Fig.[T)3-c we plot 
the dependence of the real part of the probe and coupling 
susceptibilities. 




Fig. 2. The radial profile of the probe (solid line) and of the coupling 
(dashed line) beam. We remark the different scales of the beam sizes. 
Starting at r = 200 the horizontal axis is given in logarithmic scale 
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In experiments, the spatial transverse width of the 
coupling beam is much larger that the one of the probe 
beam. Therefore, we will study the dynamics of the ini- 
tial configuration shown in Fig.O The coupling field is 
approximated by a Gaussian function of the form f(r) = 
Aq exp[— (r/w) 2 ], with maximum amplitude Aq = 100 and 
a transverse width w = 8 ■ 10 4 . Once the control field is 



properly defined, the probe beam from Fig. [5] is computed 
as a stationary state of Eq. (|la|) using a standard shoot- 
ing method, assuming a spatially constant coupling field 
Q c = 100 in the vicinity of the origin. 

3. Mapping functions 

Due to their high accuracy and facility to accommodate 
mapping functions, we choose to discretize the spatial co- 
ordinates using a Chebyshev pseudospectral method. In or- 
der to properly implement such a method, our infinite do- 
main of interest is first truncated (in each spatial direction) 
to the interval [— L, L], L = 5 • 10 5 , and then scaled (with- 
out loss of generality) to the interval [-1,1]. This scaling of 
domains allows the direct use of the Gauss-Lobatto points 
given by 



A mapping function g is defined as 
x = g(s,a), 



(4) 



(5) 



where x represents the physical coordinate, — 1 < s < 1 is 
the computational coordinate (discretized by the Gauss- 
Lobatto points), and a denotes one or possibly more free 
parameters. These new sets of collocation points s gener- 
ated through mappings of the Chebyshev points retain the 
ease of evaluation of the successive derivatives of a given 
function. For instance, the first and second derivatives of 
u(x) can be straightforwardly evaluated as 



1 



du 



du 

dx g'{s, a) ds ' 

d 2 u 1 d 2 u g"(s,a) du 

dx 2 [g'(s, a)) 2 ds 2 [g'(s, a)] 3 ds ' 



(6a) 
(6b) 



For more information related to the use of mappings func- 
tions, we refer the reader to Ref. [5]. 

The profile of our narrow probe beam (see Fig. [2]) exhibits 
an almost flat region around x = before starting its decay 
to zero. We would like to have its whole support properly 
discretized, if possible with an almost uniform distribution 
of points in order to capture all the possible dynamics that 
might take place along its spatial extent. To this intent, we 
introduce the following polynomial mapping 



x = (as + s 2 p +1 )/(1 + a), 



(7) 



where a,p > 0. Adjusting the parameters a and p one can 
control the size of the region of uniformly distributed points 
and the number of points located in this region. An almost 
uniform distributed points near the origin is achieved due 
to the nonvanishing first derivative of the mapping function 
g'(Q, a) = a/(l + a). Hence, the choice of the parameters 
a and p have to ensure that, near the origin, the dominant 
contribution comes from the first order term. Polynomial 
mappings similar to ([7]) were used in compresible mixed 
layer computation [12] in order to compare several error 
functionals of an adaptive pseudospectral method. 



Polynomial mapping, a=0.5 



-1.0 
1.0 



0.0 
-0.5 



P=1 
p=5 
p=9 
p=21 



Tan mapping 



Polynomial mapping, a=5e-3 





1.0 -1.0 -0.5 0.0 



Fig. 3. Left column: polynomial mapping J2J for different values of 
slope parameter a and polynomial order p. Right column: "tan-" and 
"sinh-" mappings iJSJ— ((9J for different values of control parameter e. 
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Fig. 4. Size of the spatial discretization step for the polynomial (x), 
"sinh-" (+), and "tan-" (o) mappings: (left) entire computational 
domain; (right) central part of comparable size with the narrow 
probe beam. The mapping parameters used to generate the grids are 
a = 5.5 ■ 10~ 4 , p = 15 for the polynomial mapping, s = 7.3 ■ 10 — 5 for 
the sinh-mapping, and e = 2.2 ■ 10~ 4 for the tan-mapping. N = 351 
in all situations. 

We will compare the polynomial mapping against two 
previously proposed families of mapping functions which 
also allow a concentration of collocation points in the center 
of the domain. These mapping functions are given by 



ir = £tan(stan 1 (l/e)), 
x = e sinh(s sinh -1 (1/e)), 



(8) 
(9) 



where e > 0. The mapping ([8]) was introduced in Ref. [4], 
and constructed in such a way so that the images of near 
step functions are almost linear. The mapping ([9]) has been 
recently proposed [27] for the study of shock waves and 
blow-up phenomena. To get more insight into the proper- 
ties of the mapping (JT])-©, we plot them and their spatial 
step size along the whole computational domain, see Fig. [3] 
and Fig. [H respectively. Optimal parameters are chosen for 
all mappings as it will be discussed below. It can be ob- 
served that both the "tan-" and "sinh-" mappings produce 
nonuniform step sizes close to x = 0, whereas the poly- 
nomial mapping is able to produce a discretization grid of 
almost constant step size in the whole central region. 



3.1. Selection of mapping parameters 

The aim of quantitatively assessing the usefulness of a 
certain mapping applied to a particular problem has been 
widely addressed in the literature [1,2,4, 11]. We follow 
here the procedure presented in Ref. [4]. Mappings 0- 
((9]) are functions of one or two parameters which are to be 
determined. As criteria we will use the functional I2 [4, 12], 
and the L2 and norms of the error 
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where L = vl — s 2 d/ds. The functional I2 represents an 
upper bound of the error made when a function is approx- 
imated using the first N terms of its Chebyshev expan- 
sion [12]. The quantity 1% offers a mapping independent cri- 
teria. The formulas (|10b[) and (|10c|) compare the N points 
polynomial interpolation of the function / against the M 
points one on a larger grid of points, i.e., N < M, hence 
being the M points polynomial interpolation taken as the 
"exact" reference. All integrals are computed using Gauss- 
Lobatto quadrature formulas. Optimal values for mapping 
parameters arc then selected in order to minimize the above 
mentioned quantities. 

Our test cases will be conducted in one dimensional 
space. Nevertheless, as our two dimensional mesh is just 
the tensor product of the one dimensional grid, the con- 
clusions from the one dimension problem can be straight- 
forwardly extended to the 2D configuration. The top-flat 
profiles found in the cubic-quintic NLS model are very 
well approximated by supergaussian functions of the form 
f(r) = Aq exp[— (r/w) 2m ] [9]. The narrow probe beam pro- 
file depicted in Fig. [5] can therefore be correctly fitted to 
this type of profiles, with fitting parameters Ao ~ 21.198, 
w ~ 1.099 • 10 2 , and m ~ 4.545. We will hence use this 
supergaussian profile as our test/input function. 

As shown in Fig. [5] for a number of discretization points 
N = 351, the quantities defined by relations (|10a[) - (|10cp 
are computed as functions of the different mapping param- 
eters. It was found that, in general, a good mapping will 
minimize both I2 and L2 quantities at the same time [4]. 
Optimal values of the mapping parameters were then cho- 
sen to minimize the L2 norm of the approximation error, 
but always comparing the shape of this functional to the 
ones of I2 and L^, hi order to ensure that these function- 
als also attain close to minima values. This choice of cri- 
teria was motivated for the unsatisfactory behavior of the 
functional I2 for the "sinh-" and "tan-" mappings for small 
values of parameter e (due to a poor discretization of the 
supergaussian profile), as well as for the infinite value of 
the derivatives of the "tan-" mapping ata; = ±lase^0 
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Fig. 5. Errors in the approximation of the supergaussian profile with 
the different mappings. Left column: errors using the polynomial 
mapping JT}. Right column: errors using the "sinh-" (solid line) and 
"tan-" (dashed line) mappings ((Sj — ([Qj . N = 351 and M = 85f in all 
situations. Note the presence of two different scales for the figures 
on the right column, the left for the "sinh-" and the other for the 
"tan-" mappings. 

(see Fig. |3|) . In addition, the functional exhibits in some 
situations a much bigger oscillatory behavior than the L2 
norm, which also makes its use more difficult for the proper 
choice of the "optimal parameters" . 

Optimal parameters for the correct discretization of the 
probe field, together with the corresponding values of cri- 
teria functions (|10ap - (|10cp . are given in Table Q] for the 
different mappings under study and for two distinct num- 
bers of discretization points, N = 121 and N = 351. The 
standard unmapped Chebyshev method is also included for 
completeness. In the case of N = 121, the functions I2, L2 
and Loo exhibit similar shapes to those shown in Fig. [5j but 
with sharper minima due to the smaller number of sample 
points. In all situations, our polynomial mapping is found 
to outperform the results obtained using the other mapping 
functions due to its ability of generating an almost uniform 
discretization grid in the whole extent of the narrow beam. 
In addition, it is noteworthy to remark that the values of 
optimal parameters a and p are noncritical. Similar results 
are obtained when compared to other mappings found in 
the literature, such as those described in [6, 18]. 

From the results presented in Table [U it can be inferred 
that the polynomial mapping is much more accurate 
than the "sinh-" mapping even when using optimal values 
for parameter e, because the latter produces much bigger 
step sizes close to the origin. Furthermore, for the "sinh-" 
mapping the I2 functional does not seem to behave as an 
upper bound of the L2 and norms, as it was previously 
demonstrated in Ref. [12]. This points out a possible poor 
discretization of the function under representation. In fact, 
the number of points has to be increased till TV = 551 in 
order to have these inequalities satisfied when using this 



Table 1 

Error comparison for the probe field when using the polynomial (0, 
"tan-" (£Sj> and "sinh-" mappings. M = 851 in all situations. U 
denotes unmapped. 

Optimal parameters 
Mapping ^ _ I 2 L 2 L x 
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NA 3.2627 20.918 
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£ 


~ 2.2320e-04 


1.6753e-04 1.4671e-10 3.6194e-10 
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~ 7.2731e-05 


2.7488e-03 4.9130e-03 1.6885e-02 


U 






NA 1.7230 18.882 



mapping. The same happens when using the "tan-" map- 
ping and a small number of discretization points (N = 
121). The value of functional 1% is not assigned (NA) for 
the unmapped Chebyshev method since in this situation 
the probe field is discretized by a single collocation point. 

However, our system of interest consists in two coupled 
beams, and therefore the coupling field has to be also 
properly discretized for our choice of mapping parame- 
ters. Table [2] presents values of functionals (|10ap - (|10c[) for 
the coupling field for the choice of parameters that best 
discretizes the narrow supergaussian profile. Even with 
a reduced number of collocation points (N = 121), the 
polynomial mapping is able to produce a fairly good de- 
scription of this field, and of comparable accuracy to the 
best of the other mappings when the spatial resolution is 
increased (N = 351). On the other hand, the "tan-" map- 
ping is not capable of describing correctly this wider pro- 
file, since it concentrates almost all discretization points in 
the center of the interval. The "sinh-" mapping, as well as 
the unmapped Chebyshev method, is able to discretize the 
control field, but was not able to represent appropriately 
the narrow probe field. 

4. Numerical simulations 

The propagation of the probe and coupling fields is sim- 
ulated using a split-step mapped pseudospectral method 
as the one presented in Ref. [21]. The linear step (Laplace 
operator) is integrated by using exponential integration of 
the transformed Chebyshev matrix, whereas the nonlinear 
step is performed by using explicit midstep Euler method. 
In order to ensure transparent boundary conditions, we 
have placed an absorbing potential to get rid of the poten- 
tially outgoing radiation [21]. Using this numerical scheme 
we have simulated the time evolution of the initial probe 
and coupling fields shown in Fig.[2j given by the NLS sys- 
tem ((T]), for all the three mappings given in the previous 
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Table 2 

Error comparison for the coupling field when using the polynomial 
J7), "tan-" JS) and "sinh-" (O mappings, using the sets of parameters 
which give optimal description of the probe field. M = 851 in all 
situations. U denotes unmapped. 

Optimal parameters 
Mapping _ . I 2 L 2 Loo 
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Fig. 6. Field amplitudes for t = 3600 computed with the polynomi- 
al-mapped Chebyshev grid, with N = 121, a = 5 ■ 10~ 4 and p = 12. 
Upper (lower) row shows the probe (coupling) field, while the left 
(right) column depicts the spatial profiles on physical (computa- 
tional) domain. 

section. The parameters of the mappings were kept hxed 
during the time evolution. The time step and the number 
of sample points are set to At = 0.1 and N — 121, respec- 
tively. As the initial fields do not constitute a stationary 
solution of the coupled NLS system (fTJ) , they will change 
their shape in the course of the numerical simulation. We 
have verified that the computational results shown bellow 
are not altered when changing the size of time step, e.g., 
At = 0.01 or 1. 

In Figs. [S][S] we plot the spatial profiles of the probe and 
coupling fields on both the physical and computational do- 
mains. Around t ~ 3600 the dynamics shows the developing 
of a peak into the coupling beam f£ c , of comparable spatial 
width with the narrow probe beam, while the probe field 




Fig. 7. Field amplitudes for t = 3600 computed with the 
"sinh"-mapped Chebyshev grid, with N = 121 and e = 7.2731 ■ 10~ 5 . 
Upper (lower) row shows the probe (coupling) field, while the left 
(right) column depicts the spatial profiles on physical (computa- 
tional) domain. 




Fig. 8. Field amplitudes for t = 600 computed with the "tan" -mapped 
Chebyshev grid with TV = 121 and e = 4.5217 ■ 10 — 4 . Upper (lower) 
row shows the probe (coupling) field, while the left (right) column 
depicts the spatial profiles on physical (computational) domain. 

only exhibits slight modifications of its spatial profile. In 
the case of the polynomial-mapped Chebyshev grid, both 
the probe and coupling fields show smooth variations in 
the associated computational domain, with their peaks and 
spatial decays correctly sampled. In especial, note how the 
almost singular structure that represents the probe field is 
perfectly approximated by this mapping even using a small 
number of grid points (JV = 121). On the other hand, the 
use of the "sinh" -mapped Chebyshev grid leads to a merely 
rectangular probe profile O p , with a poor sampling of its 
spatial decay, see the upper-right plot of the Fig. [7] This 



fact is also manifested on the peak located in the center of 
the coupling beam, see the lower-right plot of Fig. [7] 

In the case of the "tan" -mapped Chebyshev grid, see 
Fig. [HI due to its poor spatial discretization, the coupling 
beam is quickly polluted, by t ~ 600, with significant errors. 
These errors are coupled back into the probe beam which 
shows a noisy spatial profile. Hence, the subsequent time 
development of the system is altered. 



5. Conclusions 

In order to study the propagation of two coupled beams 
exhibiting spatial widths of several orders of magnitude of 
difference, we have used a two-parameter polynomial-type 
mapping function especially suitable for its use in conjunc- 
tion with Chebyshev pseudospectral methods. Using error 
criteria related to the spectral accuracy, we have compared 
the approximation error attained by the polynomial-type 
mapping against the ones obtained using previously defined 
mappings proposed to capture collapse or shock wave phe- 
nomena. We have also performed numerical simulations of 
two coupled beams propagating through an atomic coher- 
ent media, where the propagation is described by a system 
of two coupled NLS equations. While the "sinh" -mapping 
and "tan" -mappings only offer proper discretizations of the 
coupling and probe beams, respectively the polynomial- 
mapping is able to capture simultaneously all the physical 
features of both fields, still using a relatively small number 
of discretization points. The results from the comparison 
of the error criteria presented in Section [3] are also sup- 
ported by numerical simulations. Furthermore, the results 
presented in Fig. [5] indicate that the optimal values of the 
polynomial-mapping parameters are noncritical. 

It is worth emphasizing the easiness of implementation of 
the proposed mapping in comparison with the implemen- 
tation of cither a multiple scale or domain decomposition 
method. In addition, a third parameter, corresponding to 
the center of the uniform discretized region, can be easily 
accommodated into the polynomial mapping, allowing the 
tracking of moving and interacting structures of small spa- 
tial size. 
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